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I .  INTRODUCTION 


This  study  is  a  sequel  to  earlier  studies  [1-2]  performed  at  Aerojet  Electronic  Systems 
Plant  under  a  previous  contract  with  the  Office  of  Naval  Research  to  account  for  the  contribution 
of  mult^le  scattering  to  microwave  brighmess  temperature  of  sea  ice  and  its  snow  cover.  The 
application  of  the  distorted  Bom  approximation  in  conjrmction  with  the  strong  fluctuation  theory 
has  led  to  calculated  values  for  brighmess  temperatures  in  agreement  with  measurements  from 
sea  ice  and  its  snow  cover  at  frequencies  below  37  GHz  [3-4].  However,  at  frequencies  above  37 
GHz,  the  distorted  Bom  approximation  overestimates  the  brighmess  tenperature.  This  is  due  to 
the  fact  that  the  distorted  Bom  approximation  accounts  only  for  the  singly  scattered  field. 
Physical  and  mathematical  arguments  given  in  [1-2]  showed  that  inclusion  of  higher  order 
multple  scattering  would  lead  to  mote  reasonable  values  for  the  brighmess  temperature  at  higher 
frequetKies  where  die  multiple  scattering  is  expected  to  be  dominant. 

An  extension  of  the  distorted  Bom  approximation  was  carried  out  [1]  to  include  higher 
order  multiple  scattering.  Strong  flucmation  theory  was  enployed  to  derive  a  set  of  integral 
equations  describing  the  second  moment  of  the  multiply  scattered  field  for  a  multi-layered 
anisotropic  medium  with  planer  interfaces.  Various  numerical  considerations  and  the 
programming  of  these  equations  on  a  digital  computer  to  obtain  numerical  results  were  left  to  a 
further  smdy. 

The  present  study  addresses  those  numerical  problems  and  describes  results  obtained 
from  a  computer  program  written  to  calculate  brighmess  temperature  of  snow  over  ground.  To 
calculate  emissivities  based  on  the  equations  in  [1],  a  number  of  five-fold  integrals  (three  to 
calculate  second  moments  of  the  electric  field  and  two  to  calculate  the  bi-static  cross  sections 
from  these  second  moments)  needed  to  be  evaluated  in  addition  to  solving  a  coupled  set  of 
integral  equations.  This  is  clearly  a  large  numerical  problem  requiring  extensive  computer 
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resources.  Fortunately,  it  was  found  possible  to  perform  two  of  these  integrals  analytically,  thus 
reducing  the  problem  to  three-fold  integrations  and  reducing  the  computational  burden 
considerably.  In  addition,  maximum  use  was  made  of  symmetries  in  the  kernels  in  order  to 
simplify  the  integral  equations  to  be  solved. 

To  meet  the  objective  of  this  study,  Peake’s  formula  is  used  in  Sec.n  to  relate  the 
brighmess  temperatures  to  the  bi-stadc  cross  sections  associated  with  the  multiply  scattered  field. 
Then  the  synunetry  characteristics  of  the  medium  are  employed  to  simplify  the  relation  between 
the  bi-static  cross  sections  and  the  field  second  moment.  The  integral  equations  describing  the 
field  second  moment  are  stated  in  Sec.in.  The  equations  contain  kernels  with  tr^le  integrals. 
One  integral  is  over  the  medium  depth.  The  second  integral  is  over  the  specmun  of  the 
azimuthal  wave  number,  and  the  third  is  over  the  difference  between  the  heights  of  the  Green’s 
function  source  points.  Arudytic  evaluation  for  the  integral  over  the  higher  part  of  the  azimuth 
wave  number  spectrum  is  performed.  Also,  the  integral  over  the  difference  between  the  heights 
of  Green’s  fimction  source  point  is  carried  out  analytically.  The  integrals  over  the  lower  part  of 
the  wave  number  spectrum  and  over  the  medium  depth  are  carried  out  numerically. 

In  Sec.IV,  the  procedures  applied  to  solve  the  integral  equations  are  described  and  some 
limitations  of  those  procedures  are  pointed  out.  Finally  in  Sec.V,  numerical  results  for  snow 
over  a  soil  surface,  illustrating  the  behavior  of  the  brighmess  tenqjerature  as  function  of  medium 
depth  for  different  radiometer  frequencies  and  polarizations,  are  given.  Those  results  are 
compared  with  corresponding  calculations  obtained  via  the  distorted  Bom  approximation . 

In  addition  to  the  work  explicitly  contained  in  this  report,  the  contract  called  for  a  critical 
review  of  chq>ters  in  a  plaimed  monogrtq^h  on  sea  ice  (subsequently  published  in  1992  by  the 
American  Geophysical  Union  as  “  Microwave  Remote  Sensing  of  Sea  Ice”,  edited  by  F.  Carsey). 
Two  separate  reviews,  one  for  the  preliminary  draft  of  (Th^ter  10  (dealing  with  theory)  and  one 
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reviewing  the  preliminary  draft  of  the  monograph  as  a  whole,  were  written  and  provided  to  tl» 
editor  during  the  course  of  preparation  of  the  monogr^h.  A  discussion  of  this  work  is  not 
included  here. 
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n.  THE  BRIGHTNESS  TEMPERATURE  ASSCX3ATED  WITH  MULTIPLE 

SCATTERING 


The  geometiy  of  the  problem  under  consideration  is  shown  in  Fig.l,  where  the  snowpack 
is  treated  as  an  isotropic  multi-layered  medium  with  plane  interfaces  over  a  homogeneous  half 
space.  The  homogeneous  half  space  stands  for  the  ground.  Since  there  is  no  adequate 
information  about  the  physical  temperature  profile  within  the  snowpack,  the  latter  is  assumed  to 
be  isothermal.  Under  this  restriction  the  medium  emissivity  can  be  related  to  the  random  field  bi¬ 
static  cross  sections  through  Peake’s  formula  [3]; 


e.(ifco)  =  l-iRa(iko)l  )} 


(1) 


where  (  ©o,  <|>o  )  is  the  observation  direction,  a  (a  =  v,  A)  is  the  sensor  polarization,  /?«(*„) 
is  the  reflection  coefficient,  and  y^(k„ik)y  ig  =  v,  /i)  is  the  bi-stadc  cross  section  of  die  random 

field  in  the  ko  direction  with  polarization  a  due  to  a  plane  wave  propagating  in  k  (6,  ()>)  direction 
exciting  the  medium.  To  examine  the  contribution  of  the  multiple  scattering  to  the  brightness 
temperature  we  write  the  random  field  as 

Era'~~  E  n’k  E  rm 

(2) 


where  £„  is  the  random  field  due  to  single  scattering  (distorted  Bom  approximation),  and  Em 
is  the  random  field  due  to  multiple  scattering.  Since  the  bi-static  cross  sections  are  proportional 
to  the  intensity  of  the  scattered  random  field  from  (2)  we  obtain. 


4 


<  I?,  > 

Accordingly  the  bi-static  cross  section  can  be  writtm  as 
y^iko,k  )  =  Yl,(jFo.*  )  +  yi,(*o.*  ) 


(3) 


(4) 


The  first  term  in  (4)  stands  for  the  bi-static  cross  section  obtained  via  distorted  Bom 
iqjproximation  (single  scattering).  The  second  term  is  the  bi-static  cross  section  associated  witii 
the  higher  orders  of  multiple  scattering.  By  substituting  (4)  into  (1)  the  emissivity  is  reduced  to 


=  sined0d«j>{Y^(fco.fe  Hfjkork  )} 


(5) 


where , 


eiiU  =  1  - ^  sinWe  ^  yUik.,k) } 

is  the  emissivity  obtained  via  the  distorted  Bom  approximation  (single  scattering). 


(6) 


Then  by  adding  the  effect  of  the  earth  atmosphere  to  (S),  the  brighmess  temperature  of  the 
radiation  leaving  the  eardi  surface  at  z  =  0  in  direction  with  polarization  "a”  is 

TAU  =  THkJ-r  iko) 

<* 

(7) 
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where  Taiko)  brightness  tenq)erature  obtained  via  distorted  Bom  af^roximation  (single 
scattering), 


rf(jF„)=  e^^ik,)T,  +  \RAkofT^ik,) 

)  +  yi,iko^k  )]T^ik)smQdQd^ 


(8) 


and  Taiko)  ia  the  mult4>le  scattering  contribution  to  the  brightness  ten^rature, 

^"^'aSi^o^k  )  +  y’Jk,,k  )){T,-T^{h)smMQd^ 

(9) 

In  (8)  and  (9)  7,  is  the  ground  physical  temperature  and  T^iko)  is  the  brightness 
temperature  of  sky  radiation.  In  this  study  the  ground  physical  temperatures  is  taken  to  be 
independent  on  6  and  The  sky  temperature  is  assumed  to  be  a  function  of  6  only. 

As  we  can  see  from  (9),  the  contribution  of  multiple  scattering  to  the  brighmess 
temperature  depends  mainly  on  the  difference  between  the  ground  and  sky  ten^ratures.  At 
microwave  frequencies,  the  ground  physical  temperature  is  higher  than  the  sky  brighmess 
tenq)erature.  Accordingly,  the  higher  order  mult^ly  scattered  field  tends  to  reduce  the  brighmess 
temperature.  The  reduction  is  proportional  to  the  bi-static  cross  sections  associated  with  the 
multq’ly  scattered  field. 

Eq.  (7)  gives  the  brighmess  temperature  at  the  earth  surface.  When  this  radiation  is 
propagated  through  the  atmosphere,  it  undergoes  some  modifications  due  to  the  atmospheric 
attenuation  and  emission.  These  efrects  can  be  accounted  for  following  the  same  procedures 
repotted  in  [2]. 
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In  (10),  the  bi-static  cross  sections  are  related  the  field  second  nunnent  tbrouj^  the  relation  [1]; 


[it*cose/(iccose„)]  b.bj  \dz'  Idz"  Aa,iz\k^,ky) 

A*jf(z" ,k,,ky)  F^ik,,k,,z'  ,z") 


(10) 


where  symbol  *  denotes  coiiq>lex  conjugate,  die  Fha{kx,ky,z'  ,z" )  is  the  Jfcr  elonent  of  the  two 
dimensional  Fourier  transform 


FS^x*k,,T,T')=  f  dxdyCiy-r”\)txp(ik^x+ik,y)<E„i^)E*„(f')> 


C(|r’-r"|)  =  <^(r')^*(r”)> 

^(r')  =  AI,/(l  +  5AL] 

AL  =  kl{Kc-K') 

S  =  -ll{3klKc) 

(11) 

Here  Ko  is  the  quasi-static  dielectric  constant,  K'  [5]  is  random  dielectric  constant  of  the 
medium,  S  is  the  coefficient  of  the  delta  function  part  of  the  Green’s  function  and  C  (|r  -r”|)  is 
the  correlation  function. 

In  (10),  Aitiz'  ,kx,ky)  is  the  ik  element  of  the  two  dimoisional  Fourier  transform  of  the 
dyadic  Green’s  function  r(r,f’ )  with 
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(12) 


F(r,r’ ) » (2n)‘*  Jf/*,  ]dkyAiz’  ,kz,ky) «q>li[ifcx(JC - x* )  +  kyiy - y )  + 

and 

Jtx*  +  k^  +  k^  *  ^9* 

A(2’  » i4pp(z’  ,6)  pp  +  >4p»(z'  ,0)  p2  +  A^(2'  ,0)  ^  +  An.(z’  .6)  zp  +  A„(z’  ,0)  zz 

(13) 

In  (13),  p,  and  z  are  unit  vectors  along  the  cylindrical  coordinates,  and 
kz^  +  ky^  =  Q^,  p  =  ib„sin0 

(14) 

Due  to  the  medium  symmetry  around  z  axis,  (x,y),  (y,x),  and  (z,y)  componmts  of  the 
field  second  moment  dyad  vanish.  Thus  as  shown  in  [1], 

Yrt(jF<„f)  =  (fc9Cos®)V(»ccos0„)  Tdz’  T  dz"  A^(z’,0)A**t(z’’,0). 

[  sin*<|>  Fiui^  ,z"  )  +  cos*<|>  Fnai^'  .z”  ) ) 

(15) 

yLikoyk)  =  klKncos%„)  Jdz'  ]dz"  Upp(r',0) A*pp(z’,0).[cos^<|>Fn«(z',0) 

+  sin*<|>  F22<.(z'  ,z"  )  ]  +  cos<|>  [  Appiz'  ,0)  A*fn{z"  ,0)  Fi3«(z’  ,z" ) 

+Ap*(z'  ,0)  A*pp(z”  ,0)  FzUz'  ,z" ) )  +  /ip,(z'  ,0)  A*f»{.z"  ,0)  F^iaiz'  ,z")] 

(16) 
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Then  by  substituting  (15)  and  (16)  into  (9)  and  perfonning  the  integration  over  d^,  the 
ccmtribution  of  die  mult^ly  scattered  field  to  the  brightness  tenqierature  can  be  written  as 


11/2  0  0 

ri(jFJ  =  *J/(8ncose„)  J  sinOde  fdz'  jdz"  ,z'’)  + Fnaiz' ,z")]. 

[  cos*e  A^(z'  ,0)  A%(2'  *  .0)  +  App(z’  ,0)  A*pp(z’  .0)  ]  + 

+  2  Ap.(z'  ,0)  A*f.(z",Q)  Fs,a (z' ,z" ) }{r,- T^dc)) 

(17) 

In  (17),  the  Fu,  ‘s  (i  =  1,  2,  3)  are  the  diagonal  elements  of  the  p^(kx,ky,z'  ,z”)  tensor 
(Eq.  11)  when  the  exciting  field  is  polarized  in  "a”  direction  {a  =  v,h).  More  sinqilification  for 
(17)  can  be  aichieved  if  the  correlation  function  C  (|r  -r'|)  is  assumed  to  be  sharply  peaked.  In 
this  case  the  diagonal  elements  of  the  f^ikx,ky,z'  ,z'* )  tensor  reduce  to 

F./.(2’ .z”  )  =  < eii.iz'  )>W'iz' ,z'\kx,ky) 


<  e«(z’ )  >  =  <  £,v(z’  )£*i,(z' )  > 

W (z' ,z'\kx,ky)  =  ]dxdy air  -f"  |)exp(ik,x  +  ijfeyy) 

For  an  exponential  correlation  function  with 

C  (IF  -F  ’  I)  =  <  l^f  >  exp(-|  r-r\/£) 

the  Fourier  transform  W’  ( z'  ,z'' ,  kz,ky)  can  be  written  as  [1] 


(18) 


(19) 
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W{z\z'\k,,ky)-^<^(>2f-p-^[l  +  p\z'-z''\]exp{-p\z'-2''\^ 
where 


(20) 


(21) 


Substituting  (20)  into  (18)  then  substituting  the  resultant  in  (17)  with  and  letting  z' s  zi  and 
z"  =  zt  +  u  (  see  Appendix  2  in  [3]  for  a  discussion  of  this  transformation),  the  multiple 
scattering  contribution  to  the  brightness  temperature  can  be  written  as 

*/i  0 

ri(jtJ  =  A:2/(2cos6„)  j  sinOdO  idziRe{[£iia(zi)  + £22^(21)) [cos^0|A4»(zi)|  /#(zi,0) 

0  -« 

+|App(zi,e)  f  /pp(zi  ,0)]  +  2  £33.(zi)|Ap*(z,.0)|"  /pp(z,  ,0)  ){r,  -  T^tyik)) 

(22) 


where  Re  is  die  real  part  operator,  and 

/♦♦(zi,0)  =  ^  j  dMir(|«|,fc„ity)exp(~  |a-(z”,0)dz"  ) 

n+ii 


/pp(zi,0)  =  ^  1  duW'Hk.,ky)t^i-  ?P_(z”,0)dz'’ ) 

—  TI+« 

(23) 

Here  the  functions  ol  and  are  related  to  the  (jreen’s  function  and  were  discussed  in 
Appendix  A  of  [1].  For  convenient  reference,  they  are  reviewed  in  Appendix  A  of  this  study. 
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Explicit  expressions  for  ol  and  are  derived  in  Af^ndix  B.  The  analytic  evaluation  for  the 
integrals  in  (23)  is  performed  in  Ai^)endix  C. 

Eqs.  (22)  and  (23)  give  the  multiply  scattered  field  contritmtion  to  the  brighmess 
ten^rature  for  an  isothermal  snowpack.  As  we  can  see  from  (23)  the  brighmess  temperature 
over  the  snowpadr  surface  is  r'^ated  to  the  field  second  moments  within  the  snowpack.  In  the 
following  section  the  integral  equations  describing  second  moments  of  die  multqile  scattered 
field  will  be  considered  and  its  reduction  to  a  form  suitable  for  digital  conqiuter  programming 
will  be  investigated. 
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m.  INTEGRAL  EQUATIONS  FOR  THE  MULTIPLY  SCATTERED  FIELD  SECOND 

MOMENT 

In  [1]  two  uncoupled  sets  of  integral  equations  have  been  derived  for  die  multiply 
scattered  field  secoml  moments.  One  set  describes  the  diagonal  components  of  the  field  second 
moment  tensor,  and  the  other  set  describes  the  off  diagtmal  terms.  Since  the  brighmess 
temperature  for  a  linearly  polarized  signal  depends  only  on  the  diagonal  continents  (see  Eq. 
17),  they  are  the  only  conqionents  that  will  be  considered  here.  From  [1],  we  can  write  the 
integral  equation  describing  those  components  as 

0 

Eii«(z)  —  \dz\  {Dii(2,ztX^ii«(zi))  +  Di2(z»ziXe22*(zi))  +  bfti2,z\){e^u^{z^)) 

0 

Enti^)  =  jdzi  {Dij(2,ZiX«U<i(Zi))  +  Dll(Z,ZlX^22«(Zl))  +  b(»(2,Zl)(e3J<.(zl)) 

0 

Etia(,z)  =  jdzi  {b^(Z,ZlX«tu(zi))  +  btp(2,ZlXe22«(zi))  +  2b„(2.ZtXe33«(Zl)) 


(eiiaizi))  *  Eua(.Z\)  +  Euaizi),  1  =  1,2,3 


(24) 


where  E‘„.  are  the  second  moments  of  the  coherent  field  when  the  exciting  field  is  polarized  in 
"a”  direction,  and 


Dn(2,zi)  =  3  /  4(bpp(z,2i)  +  b^{z,zi)  + 1  /  3dpt(2,zi)  ] 


Di2(2,Zi)  =  1  /  4{bpp(z,Zl)  +  ^♦♦(2,Zi)  -  bpt(2,Zi)  ) 


(25) 
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Here 


M  0 

b^{z,zi)  =  ^  Re(  f  pdp  iduW  (|4zi)  a^{z,zi)  a*«tp(z,zi  +  u) } 

0 

(26) 


where  the  subscripts  (afi)  are  one  of  the  pairs  (pp),  (44X  (pz)>  or  (zp).  In  (25), 

m  0 

hp4(z,zi)  =  ^Re{  Jp<ip  fduW'(^^zi)[app(z,zt)  a*^(z,zi  +  u)  +  a^(z,zi)  a*pp(z,zi  +  «)]} 

0  -M 

(27) 

In  (26)  and  (27)  Uos’s  are  the  components  of  the  Green’s  function  tensor  with  the 
snowpack  which  are  given  in  ^^ppendix  A.  Their  relation  to  the  Green’s  function  elements  in  the 
upper  half  space  (Eq.  13)  is  also  given  in  Appendix  A. 

The  problon  of  integrating  over  large  values  for  the  spectrum  of  azimuth  wave  number  p 
can  be  overcome  by  writing  the  integral  over  p  as 

Jdp=Jdp  +  Jdp 

0  0  ft 

(28) 

where  Pi  is  a  very  large  wave  number.  Then,  by  making  use  of  the  Green’s  function  behavior  for 
very  large  p,  the  second  part  of  integration  over  p  can  be  carried  out  analytically  as  shown  in 
Ai^pendix  C  of  [1]  reducing  (24)  to 

0 

Eiu(z)  +  mi(z)(eiia(z))  +  m2(z)(e22a(z))  =  jdzt{  Mn(z,zi)  (cii«(zi))  +  Mi2(z,zi)  (e22«(zi)) 

-•ee 

+  /(»(^»Zl)  (e33«(zi)) } 
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0 

miiz){eiu(z))  +  E^(z)  +  miiz){en^(2))  =  fdzti  MuiZfZi)  (eij«(zi))  +  Afii(z,2i)  {ej2«(zi)) 

+  /pr(2.Zl)(<f3i-(Zl))) 


0 

+  W3(2)  (e33-(o)  =  \dz\{  f^iz,z\)  {ewAzi))  +  f^iz.zx)  {enaizi)) 

+  Afj3(2,Z,)(e33,(2l))  } 


(29) 


Explicit  expressions  for  wi(z),  mi(z) ,  and  m3(z)  are  given  in  [1].  In  addition 
Af  ii(z,2t)  =  3  /4{/pp(z,zi)  +  /4*(z,zi)  - 1  /  3/p^(z,zi) } 

Af2a(z,zi)  =  1  /  4{/pp(z,zi)  +  -  f(^(z,zi) } 

A/33(Z,Zi)  =  2  f„(Z,Zl) 

(30) 

Here  the  subserved  functions  /  are  idmtical  to  the  subscr^ted  functions  b  defined  in 
(26)  and  (27)  accept  for  the  fact  that  the  upper  limit  of  the  integration  over  p  is  bounded  by  pi 
and  the  singularity  is  extracted  from  as- 

A  further  simplification  for  the  integral  equations  can  be  achieved  if  the  integration  of  the 
subscripted  functions  over  u  is  evaluated.  This  can  be  achieved  by  using  the  translational 
property  of  the  Green’s  function  tensor  elements  (  Appendix  A)  leading  to; 

/4^(z,zt)  =  ^Re{  fpdp  ^duW i\u\,zi)\a^iz,zi)\^3*^{zuu)] 

0 
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/pp(z,2i) »  ^  Re(  Jpdp  fdu  W  (|4z,)  |  a^(z,zi)  f  3*pp(zi,«) ) 

0 

/p,(2,zi)  =  ^  Re{  /p</p  /</«  W (14z,)  I  a^{z,zi)  | ^^*„{zuu) ) 

0 

/,p(z.zi)  =  in  Re{  Jpdp  \du  W 0421) I  a^{z,zi)  |^3*pp(2i,M) ) 

0 

/„(2,ri)  =  ^  Re{  ipdp  fdu  W (|u|,2i)  1  «„(2,zi) \^S*„izi,u)  ] 

0 

/p4(z.2i)  -is  Re  /dp  jduVr  (|tt|,zi)  { app(2,zi)  <i*h(z.2i)  3*h(z.2i) 

0 

+  a#(2,Zl)a*pp(2,Zi)  3*pp(2,Zl)  ) 


where 


5^(zuu)  =  exp[-  fa.(z’ ’)dz”] 


n*M 

«+■ 


=  exp[/ou(2”)d2”] 


XI 


and 


3pp(2i,ii)  =  exp[- jf  P.(2")d2''] 


xi+» 

*!+■ 

=  exp(  /P^(2”)d2”] 

XI 


if  r  >  zi  +  « 
if  2<Zl  +  « 

if  2  >  Zi  +  M 
if  z<zi  +  u 


Then  by  carrying  out  the  integral  over  u  in  (31)  we  obtain 


(31) 


(32) 


(33) 
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Pi  2 

/»(i,i.)=  lpdp\  Rc{/4«(zi,P)) 

0 

/pp(z,zi) *  /  p</p|  app(z,zi)  I*  Re{  /pp(zi,p) } 

0 

/„(z,zi)  =  J  pdpl  ap,(z,zi)  p  Re{  ] 

0 

/,p(z.zi)  =  J  p<ip|  a^(z,zi)  I*  Re(  /pp(zi,p) ) 

0 

/„(z,zi)=J  p</p|  a„(z,zi)  p  Re{/pp(zi,p)} 

0 

Pi 

/p*(Z.Zi)  =  /  prfpRe  {  dpp(z,Zj)<I*44(Z,Zi)  I*^(zup)+  ^♦♦(z,Zi)a*pp(Z,Zi)  /*pp(zi,p) 

0 

(34) 

where 


/h(zwP)  =  ^  (|«Ui)  3*t(zi,p) 

»«• 


/pp(zi,p)  =  ^  J  W’  (Hzi)  3pp(zi,p) 


(35) 


An  implicit  dependence  of  and  /pp  on  z  has  been  suppressed  in  (35). 


For  the  exponential  correlation  function  given  in  (19),  an  explicit  expression  for  can 


be  obtained  from  Appendix  C  as. 


<14f  > 


Zc{p,a)  +  z,ip,a)  I 


_<l^r>  f  - _x  .  a+(zi)  _  _x  .  [a+(zi)-a+(zi)] 


where 


Z  >  Z) 


z<  z, 
(36) 


Zc(p,a)* 


2 

[p*  +  a(zi)]* 


Zsip,a)  =  -^ja{zl) 


3/>^  +  a(zi)] 
p*(p^  +  a(2i)] 


3p^  +a(zi) 

p^  +  aizi) 


+  p{z 


i~z)  • 


sm[V<i(zi)(zi  -  z)] + cos[Va(zi)(zi  -  z)3  -  ^  sin[Va(zi)(zi  -  z)] ) 


(37) 


The  functions  a(zi)  are  related  to  Green’s  function  tensor  and  are  defined  in  Appendix  A. 
Similar  expressions  can  be  written  for  /pp  . 


Eqs.  (29)  are  the  reduced  form  of  the  integral  equations  describing  the  second  moment  of 


the  multiply  scattered  field.  In  the  next  section  they  will  be  put  in  a  form  suitable  for  solution  by 
digital  con^uter. 


IV.  SOLUTION  OF  THE  INTEGRAL  EQUATIONS 


To  put  Eqs.  (29)  into  a  fonn  suitable  for  digital  computer  programming,  we  transfer  than 
to  a  set  of  linear  algebraic  equations  which  can  be  arranged  into  a  linear  matrix  equation.  Then 
matrix  inversitm  techniques  can  be  qjplied  to  solve  the  resultant  matrix  equation. 

To  transfer  Eqs.  (29)  into  a  set  of  linear  algebraic  equations,  the  integration  over  p  and 
over  zi  must  be  evaluated.  Gaussian  quadratures  [6]  will  be  enq>loyed  here  to  evaluate  the 
integrals  over  p. 

As  for  integrals  over  Zi  we  let 

0  0 

i  !  dz\+  /  dzi 

-g, 

(38) 

where  z$  is  the  snow  depth.  Since  the  ground  beneath  the  snow  is  homogoieous,  die  integration 
kernels  vanish  for  any  source  point  located  below  z,.  Accordingly,  the  second  part  of  integration 
in  (38)  reduces  to  zero,  and  the  lower  limits  of  integration  over  zi  in  Eqs.  (29)  can  be  r^laced  by 
z  =  z,  ■ 

The  snow  dq>th  z,  is  divided  into  (N-1)  equal  increments  of  dz  each.  Then  by  qrplying 
the  trapezoidal  rule,  the  integration  over  dz  in  (29)  turns  into  a  summation.  The  trapezoidal  rule 
is  tq^lied  here  because  through  this  rule  the  source  point  of  the  Green’s  function  can  be  located 
(z  =  Zi)  and  the  discontinuity  in  Sp,  and  a^,  at  this  point  can  be  handled. 

By  introducing  the  vectors 
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£ii«(z.) 

EtlmiZm) 

Ess^iz,) 


£fi.(z,) 

£52.(2.) 

£52.(2.) 


Eqs.  (29)  can  be  written  as 

7^«  +  W4£«  +  ^~]=i  Wnj[Ej,  +  E%] 


where  /  is  the  unit  matrix  and 


(39) 


(40) 


m, 


mi  (2.)  m2(z.)  0 

m2 (2.)  mi  (2.)  0 

0  0  m3(2«) 


Afll(2«,2y)  Mi2(z.,Z;)  /p,(2«,2;) 

Mi2(2«,2;)  Mii(z.,Z;)  /p,(z.,Z;) 

/^(z  .♦2y)  ^  2p(Zff9Zy)  ^jy(z.jZy) 


(41) 


Eq.  (40)  represents  a  set  of  3N  algebraic  equations.  To  solve  these  3N  equations  we 
introduce  new  column  matrices  £‘  having  3N  components.  The  elements  of  these  two 

new  vectors  are  related  to  the  elements  of  the  field  vectors  given  in  Eq.  (39)  through  the  relations 
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£-(0*  £-(3(n-l]  +  i), 


£5.(0  =  K(3[n-l]  +  O. 


i*  1.2,3 


Accordingly,  the  3N  algebraic  equations  resulting  hxnn  (40)  can  be  arranged  as  follows: 


[7-a]e,=a£: 


where  A  is  a  square  matrix  of  order  3N.  Its  elemoits  are  related  to  the  elements  of  the  matrix 
givoi  in  (42)  through  the  following  relation : 


mn,(i,/)5»»  +  Ar*,((,y)*A({3r/t-l]+/},(3(/«-l]+y7),  i,  j  =1,2,  3 


Eq.  (43)  is  the  matrix  equation  governing  the  second  moment  of  the  multiply  scattered 
field.  It  may  be  solved  by  using  either  eigen  value  methods  or  a  matrix  inversion  technique. 
Since  the  eigen  values  technique  requires  more  computer  resources  (memory  and  computing 
time),  only  the  matrix  inversion  techniques  was  considered  in  diis  study. 

According  to  matrix  inversion  technique  we  premultiply  Eq.  (43)  from  both  sides  by  the 
inverse  of  the  matrix  [7  -  A]  to  write  the  multiply  scattered  field  as 


£.  =  [7- A  ]"ab'. 
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The  solution  given  in  Eq.  (45)  is  valid  only  wh«i  the  determinant  of  the  matrix  [/  -  A] 
has  values  other  than  zero.  Numerical  calculaticms  using  the  developed  FORTRAN  program 
showed  that  die  determinant  of  the  matrix  [/  -  A  ]  is  controlled  by  two  factors: 

(i)  the  operating  frequency  ( operating  wavenumber  k«),  and 

(ii)  die  snow  depth. 

The  determinant  has  a  maximum  value  equal  to  unity.  This  value  occurs  when  either  the 
operating  frequency  or  the  snow  depth  is  equal  to  zero.  In  both  cases,  there  is  no  scattering 
occurring  within  the  medium,  and  the  elements  of  the  matrix  A  reduce  to  zero.  Increasing  eidier 
the  operating  frequency  or  the  snow  depth  leads  to  an  increase  of  tlie  multqile  scattering 
(increasing  elements  of  A ).  The  latter  reduces  the  determinant  until  it  reaches  zero.  In  this  case 
the  operating  frequency  reaches  a  value  equal  to  one  of  the  iMdium  resonance  frequencies  (eigen 
values).  For  frequencies  higher  than  that  frequwicy,  Eq.  (45)  gives  norqihysical  values  for  the 
multiply  scattered  field.  It  may  lead  to  a  scattered  intensity  higher  than  the  incident  intensity  or 
may  give  an  intensity  having  negative  values. 

Based  on  the  above  discussion,  the  solution  in  (45)  for  the  multiply  scattered  field  is  valid 
only  when  the  determinant  of  the  matrix  [/  -  A]  has  values  greater  than  zero.  This  can  be 
achieved  by  lowering  either  the  operating  frequency  or  decreasing  the  snow  depth.  The 
following  table  shows  the  maximum  snow  depth  to  which  solution  in  Eq.(45)  can  be  ^lied 

Frequency  (GHz)  37.0  50.0  85.0  91.6 

snow  dq>th(cm)  12.0  7.50  3.45  2.50 


Table  I 
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Table  I  is  calculated  for  uniform  snow,  with  snow  density  equal  to  0.2S  g/cm’  and  grain 
diameter  equal  to  0.7mm.  Numerical  calculations  show  that  the  values  of  the  snow  depth  given 
in  Table  I  can  be  increased  by  increasing  the  snow  density  or  reducing  the  upper  limit  of 
integration  over  p  (pi  in  Eq.  (29) ).  For  values  of  Pi  greater  than  ko ,  total  reflection  occurs 
within  the  snow  layer  leading  to  a  surface  wave  propagating  along  the  layer  interfaces. 

The  numerical  limitations  in  solving  the  equation  for  the  field  second  moments  in  the 
form  discussed  above  were  discovered  rather  late  in  the  course  of  this  work.  Overcoming  these 
restrictions  will  require  significant  modification  in  the  FORTRAN  program.  Thus  the  examples 
to  be  discussed  in  Sec.  V  of  the  report  will  not  cover  the  full  range  of  conditions  of  interest  in 
studying  naturally  occurring  snowpacks. 


The  following  paragraph  will  indicate  some  measures  that  can  potentially  be  useful  in 
extending  the  range  of  validity  of  the  numerical  calculations  in  future  studies. 

One  method  to  increase  the  range  of  conditions  for  which  a  numerical  solution  is  possible 
is  to  rearrange  Eq.  (40)  as 


+  m*,)£,a  +  m»,£j 


(46) 


By  premultiplying  (46)  by  the  inverse  of  the  matrix  (7  +  m«,),  the  equation 

_  S  —1  —  **1  ^  -  -  r  __  T 

£iio  +  (  7  +  mwf  )  +mw,)  ^  M^\  Eja"^  Eja\ 

/=! 


(47) 
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is  found.  Then  using  Eq.  (47)  we  obtain  a  matrix  equation  similar  to  Eq.  (43).  The  determinant 
of  the  resultant  matrix  is  greater  than  the  determinant  of  the  matrix  [7  -  A  ]  of  Eq.  (43). 


A  further  improvement  for  the  solution  can  be  achieved  if  in  Eq.  (40)  we  isolate  the 
Green’s  function  at  the  source  point  (n  =  j)  in  one  side  of  the  equation; 


U+m„-  MjiE„  +  £L)  =  £«  +  i  M  J  £>.  +  E%] 

jmn 


Eq.  (48)  may  be  rewritten  as 


[j  5jn-il  +  M„j\(Eja'^  E‘ja)-il + 


(48) 


J0n 


(49) 


Eq.  (49)  is  expected  to  lead  to  more  stable  matrix  equation  than  Eq.  (43).  This  is  because 
the  matrix  on  the  right  hand  side  of  (49)  reduces  to  the  unit  matrix  for  lower  or  higher  values  of 
the  operating  frequency  or  snow  depth.  The  limit  at  lower  values  arises  from  the  zero  value  of 
all  mrn^tl  elements.  The  limit  at  higher  values  of  the  operating  frequency  or  the  snow 
depth  arises  because,  the  contribution  of  the  source  point  in  (49)  increases  rapidly  with 
increasing  the  operating  frequency.  On  the  other  hand  all  other  terms  decrease. 

Accordingly  the  second  part  of  the  matrix  in  the  right  hand  side  vanishes  leading  to  a  unit  matrix 
to  be  inverted.  Since  programming  Eq.  (49)  requires  considerably  more  effort,  it  is 
recommended  for  another  study. 
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V.  NUMERICAL  RESULTS  AND  DISCUSSION 


The  equations  discussed  in  the  previous  sections  were  used  to  calculate  the  brightness 
temperature  of  snow  at  various  frequencies  and  angles  for  particular  snow  conditions.  These 
may  be  conq>ared  with  the  results  obtained  using  the  distorted  Bom  q>proximation  and  reported 
in  [3].  Reference  [3]  also  shows  ranges  of  measured  brighmess  temperatures  at  various 
frequencies  at  a  50o  angle  of  incidence  obtained  from  satellite  borne  radiometers.  From  the 
discussion  in  [3],  a  reasonable  snow  model  describing  the  conditions  at  the  satellite  ground  tmth 
sites  is  snow  at  a  tenq)erature  of  -S^  C(268.15K),  density  p  =  0.2Sg/  cm^,  and  a  grain  diameter  d 
=  0.7  nun.  These  values  are  adopted  for  most  of  the  conq>utations  to  be  reported  here  together 
with  the  same  atmospheric  model  and  underlying  soil  conditions  described  in  [3]. 

Figs.  2-4  show  comparisons  of  the  new  calculations  with  the  distorted  Bom 
approximation  as  a  function  of  the  angle  of  incidence,  all  for  a  radiometer  at  satellite  altitude. 
The  curve  identifiers  are  the  same  in  all  of  the  figures.  The  initial  letter  (h  or  v)  identifies  the 
polarization  (horizontal  or  vertical)  while  the  second  letter  (m  or  s)  identifies  the  theory  : 
mspresent  theory  with  multiple  scattering,  s  =  distorted  Bom  ^jproximation  (single  scattering). 

Figs.2-4  show  the  substantial  effect  of  including  the  multiple  scattering  terms  in  the 
brightness  temperature  computations  at  high  frequencies.  These  effects  increases  as  the 
frequency  increases  from  37  GHz  (Fig.2)  to  50  GHz  (Fig.3)  to  91.6  GHz  (Fig.4)  as  one  would 
expect.  Thus,  for  example,  at  an  angle  of  incidence  of  50^  and  horizontal  polarization,  the 
inclusion  of  multiple  scattering  reduces  the  computed  brighmess  temperature  compared  to  the 
distorted  Bom  approximation  theory  by  5.4K  at  37  GHz  to  12.  IK  at  50  GHz  to  87.7K  at  91.6 
GHz.  The  corresponding  decreases  for  vertical  polarization  are  4.2K  at  37  GHz,  9.5K  at  50  GHz, 
and63.4Kat91.6GHz. 
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Another  perspective  on  these  results  is  obtained  by  examining  emissivities.  Again 
choosing  an  observation  angle  of  50*^ ,  die  emissivities  conespcmding  to  Figs.  2  -  4  are  shown  in 
Table  n.  It  is  clear  that  the  distorted  Bom  tqiproximation  emissivities  readi  a  minimum  between 
37  and  91GHz  and  then  increase  with  increasing  frequency.  This  characteristic  behavior  of  the 
distorted  Bom  {Approximation  calculation  is  not  observed  in  measured  data  and,  in  fact  was  the 
princ4>al  reason  for  looking  at  more  conqilete  theories  such  as  the  subject  of  this  report.  On  the 
other  hand,  with  multiple  scattering  included,  the  emissivities  decreases  with  increasing 
frequency. 

TABLE n 

Calculated  emissivities  of  model  snow  padk  at  6  =  50^  (p  =  0.25g/cm^,  d=0.7  mm) 


Frequency 

(GHz) 

Horizontal  Polarization 

Vertical  Polarization 

Distocled 

Multiple 

Distmted 

Multii^ 

Bom 

Scatter 

Bom 

Scatter 

37.0 

0.8509 

0.8439 

0.9126 

0.9054 

50.0 

0.8271 

0.7725 

0.8941 

0.8399 

91.6 

0.8715 

0.5352 

0.9342 

0.7315 

Fig.S  shows  computed  brightness  temperature  as  a  function  of  snow  depth  at  an  angle  of 
incidence  of  50®  for  a  frequency  of  37  GHz.  As  discussed  in  Sec.  IV,  the  present  calculation 
method  for  multiple  scattering  is  restricted  to  depths  that  are  not  too  large.  For  the  exanqile  in 
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Fig.5,  the  maximum  snow  depth  at  which  the  mult^le  scattering  results  should  be  ctmsidered 
trustworthy  is  about  12  cm.  With  this  restriction  in  mind,  it  is  interesting  to  cmnpare  the 
calculations  in  Fig.S  with  satellite  measurements.  As  discussed  in  [3],  the  reported  range  of 
brighmess  temperature  values  over  the  ground  truth  sites  was  202  to  224K  for  horizontal 
polarization  and  228  to  2S0K  for  vertical  polarization  for  snow  depths  between  5  and  20  cm.  The 
distorted  Bom  approximation  for  horizontal  polarization  yields  brightness  temperature  ranging 
from  234.4K  to  22S.4K  in  the  snow  depth  range  of  6  to  IS  cm.  This  is  on  the  high  side.  On  die 
other  hand,  with  multiple  scattering  included,  computed  brighmess  temperatures  range  from 
229.1  to  208.4K  in  the  6  - 12  cm  depth  range  (recall  the  eadier  warning  about  depths  larger  than 
qiproximately  12  cm)  and  reach  224K  at  a  depth  of  7.5  cm.  These  values  are  seen  to  fall  very 
nicely  within  the  measured  range.  For  vertical  polarization,  the  distorted  Bom  approximation 
calculations  yield  brightness  temperature  values  ranging  from  248.2K  to  240.SK  in  the  S  to  IS 
cm  snow  depth  range.  While  tending  to  be  high,  the  calculations  fall  within  the  measured  range. 
However,  the  new  multiple  scattering  theory  yields  brighmess  temperatures  in  the  244. IK  to 
227.3K  range  for  vertical  polarization  in  S  to  12  cm  snow  depth  interval.  This  is  comfortably 
within  the  measured  range. 

In  [3],  some  discussion  was  devoted  to  the  effects  of  non-uniformities  such  as  a  surface 
crust  in  the  snowpack.  At  37  GHz,  non-uniformities  in  the  snow  structure  were  necessary  in 
order  for  calculations  with  distorted  Bom  approximation  to  match  some  of  the  low  brighmess 
temperatures  that  were  measured  for  horizontal  polarization.  In  particular,  a  model  snowpadc 
with  a  l.S  mm  surface  crast  with  a  density  of  0.5  g/  cm’  and  a  grain  diameter  d  =  1  mm  was 
assumed  to  cover  an  otherwise  uniform  snowpack  with  the  parameters  used  above.  When  the 
multiple  scattering  theory  is  qrplied  to  these  conditions,  the  brighmess  temperature  for  a  snow 
depth  of  10  cm  was  found  to  be  19SK  for  horizontal  polarization.  Extrapolating  possible  effects 
to  larger  snow  depths  (recall  the  present  limit  to  about  12  cm  so  that  a  direct  computation  is  not 
possible)  indicates  that  the  calculated  brighmess  temperature  in  the  range  170K  to  180K  may  be 
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found  at  horizontal  polarization  at  37  GHz.  This  is  near  the  lower  bound  of  measurements  made 
over  Siberia  for  large  snow  depths  and  is  much  better  than  the  lowest  value  of  202K  found  with 
distorted  Bom  q)proximation  under  the  same  conditicnis.  If  lower  bulk  snow  (tensities  were  to 
be  used  together  with  a  very  cold  and  dry  atmosphere,  the  distorted  Bom  approximation  does 
lead  to  brightness  tenq)erature  with  horizontal  polarization  that  are  closer  to  the  lower  limit  of  the 
measurements  over  Siberia  as  was  shown  in  [3].  However,  the  calculated  vertically  polarized 
brightness  temperatures  are  still  t<x>  high  (lowest  value  achieved  is  214K  with  the  reduced  snow 
density  and  238K  with  p  s  O.S  g/cro?  and  the  standard  atmospheric  conditions)  conq>ared  to  the 
measured  low  of  196K.  With  the  multiple  scatter  model,  even  using  a  snow  (tensity  of  0.25  g/cm’ 
yields  a  vertically  polarized  brightness  temperature  of  208K  for  a  snow  depth  of  10  cm. 
Extrqx>lating  to  larger  snow  depths  and  lower  snow  densities  indicates  that  values  near  196K  can 
be  reached  when  multiple  scattering  is  included. 

New  information  at  even  higher  fretpiencies  has  become  available  since  the  publication 
of  [3].  This  is  the  data  obtained  by  SSM/I  instrument  at  85.5  GHz  at  an  angle  of  53.1^. 
Brighmess  temperature  measurements  over  snow  by  the  SSM/I  at  this  frequency  are  often  in  die 
150  -  180K  range  for  horizontal  polarization  although,  depending  on  snow  conditions,  values  in 
the  200  to  240K  range  are  also  found.  For  vertical  polarization,  brighmess  temperatures  as  low  as 
160  K  are  sometimes  found,  but  more  typically  are  in  die  230  to  255K  range.  According  to  Fig.4 
(at  the  slightly  displaced  fretpiency  of  91 .6  GHz),  the  calculated  horizontally  polarized  brighmess 
temperature  at  an  angle  of  53^  is  162K  when  multiple  scattering  is  included  versus  244K  for  the 
distorted  Bom  approximation  if  the  snow  density  is  0.25  g/cir?  .  Thus,  the  multiple  scattering 
model  result  fits  comfortably  within  the  measured  range  while  the  distorted  Bom  approximation 
value  is  much  too  high.  Computations  were  also  performed  for  snow  densities  of  0.3, 0.4  and  C.5 
g/cm’  at  91.6  GHz.  For  horizontal  polarization,  calculated  brighmess  temperatures  at  an  angle  of 
53®  and  using  the  multiple  scattering  theory  increased  to  177K  at  0.3g/  cm^  to  209K  at  0.4  g/cm^ 
to  224K  at  0.5  g/  cm^  showing  that  a  reasonable  range  of  snow  densities  can  account  for  the 


27 


range  of  brightness  temperatures  measured  by  SSM/I.  In  contrast,  the  distorted  Bom 
approximatitm  yields  ImrizontaUy  polarized  brightness  ten^ratures  in  the  small  range  236  - 
243K  for  the  same  snow  density  range  and  actually  decrease  with  the  increasing  snow  dmsity. 
These  are  too  high.  For  vertical  polarization,  multiple  scattering  theory  yields  brighmess 
temperatures  of  196,  206,  231,  and  246K  for  mow  densities  of  0.25,  0.3,  0.4  and  0.5  g/cm^ 
respectively.  This  is  nicely  within  the  range  of  values  measured  by  the  SSM/I.  On  the  other 
hand,  the  calculated  distorted  Bom  approximation  brighmess  temperatures  are  in  the  veiy  small 
range  255  -  257K  for  vertical  polarization  as  the  snow  density  range  from  0.25  g/cm^  to  0.5 
g/cm^.  Clearly,  the  distorted  Bom  approximation  calculations  can  not  explain  the 
measurements. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  comiHjtations  discussed  in  Sec.  V  show  that  a  very  great  iminovement  in  agreement 
between  the  calculated  and  measured  snow  brighmess  temperatures  is  achieved  when  using  the 
extension  of  die  distorted  Bom  a{^roximation  (multiple  scattering)  theory  as  compared  to  the 
distorted  Bom  t^roximation  itself  at  frequencies  above  30  GHz.  This  agr^mnit  is  obtained 
without  the  arbitrary,  ad-hoc  parameter  adjustments  (often  euphemistically  called  effective 
parameters)  that  are  made  in  radiative  transfer  models  and  relies  strictly  on  random  media 
parameters  that  can  be  measured. 

It  is  particularly  gratifying  that  calculations  at  91.6  GHz  (near  the  SSM/I  8S.S  GHz 
frequency  and  coinciding  with  the  SSMA'-2  and  SSMIS  91.6  GHz  channels)  using  reasonable 
snow  parameters  yielded  brighmess  tenqieratuxes  within  ranges  measured  by  the  SSM/I  for  both 
horizontal  and  vertical  polarizations.  These  values  were  much  lower  than  those  obtained  by  use 
of  the  distorted  Bom  approximation.  The  latter  were  unreasonably  high  at  91.6  GHz. 

(Certain  numerical  difficulties  (discussed  in  Sec.IV)  were  discovered  which  prevented  the 
confutation  of  brighmess  temperatures  with  multiple  scattering  theory  over  the  full  range  of 
snow  depths,  densities  and  radiometer  frequencies  of  practical  interest.  Contract  funds  for 
exploring  possible  corrective  actions  were  depleted  when  the  problem  came  to  light.  However, 
possible  strategies  for  overcorrung  this  computational  problem  were  discussed,  and  it  is 
recorrunended  that  a  follow  on  smdy  be  undertaken  to  resolve  the  problems.  If  these  problems  are 
resolved,  the  theory  may  then  be  applied  to  the  smdy  of  multi  -  year  sea  ice  with  and  without 
snow  cover.  On  the  basis  of  the  results  obtained  to  date,  one  may  expect  to  obtain  the  first 
quantitative  explanations  of  the  high  frequency  measurements  obtained  by  the  SSM/I  over 


29 


multi  -  year  ice  as  well  as  measuiemeiits  to  be  obtained  by  future  instruments  such  as  the  SSMIS 
at  even  higher  frequencies  (e.g.,  ISO  GHz). 
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ko 


Homogeneous  ground 


Fig.l  Geometry  for  brightness  temperature  computation. 
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Fig.  2  Computed  brightness  temperatures  of  uniform  dry  snow  at  37  GHz  and  depth  of  6  cm 
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APPENDIX  A 


THE  DYADIC  GREEN’S  FUNCTION 


As  seen  in  Sec.n  and  Sec  Jn,  the  dyadic  Green’s  function  elements  aze  required  to 
determine  the  random  field  second  moment.  Since  the  medium  horizontal  dimensions  extend  to 
infinity,  it  is  convenimt  to  express  the  dyadic  Green’s  function  tensor  in  terms  of  a  two 
dimensional  Fourier  transform. 


r(r,f')  =  ^  ]dk,  ]dkyaiz,z\kx,ky)exp[j[kxix-j^)  +^y(y-y)]) 


where 


d(r,z',/t„^y)*appPP  +aprpi  +  +  a.pzp  +  aaH 


(A-1) 


(A-2) 


Explicit  formulations  for  the  Green’s  function  tensor  elemmts  in  terms  of  the  medium 
dielectric  properties  have  been  derived  in  [  1  ].  In  this  Appendix  only  the  formulations  required 
for  the  present  study  will  be  summarized. 


The  first  element  drat  will  be  considered  is  a^  which  has  the  following  formulation  [1]; 


a4*  =  exp[  Ja+(z”)d2"  ]/[a-(z')-a+(z')]  ifz>z’ 

=  exp[  -  Ja_(r'’  )dz"  ]  /  [  a-(z’ )  -  a+(z’ )  ]  if  z  <  z' 

(A-3) 
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The  a ‘s  in  (A-3)  are  governed  by  the  differential  equation 


da± /<*  +  oc|  +  a(z)  =  0 

with 

a{z)  =  klKoiz).p‘^,  =  + 

Eq.  (A-4)  is  subjected  to  the  boundary  conditions 

a+(0)  =  j 

a-  (-<»»)  =  -  iV  <*(-**) 

At  a  point  of  dielectric  constant  discontinuity  (z  =  -  d ),  a  ‘s  satisfy 


(A^) 


(A-5) 


(A-6) 


«*(-</)  c:;=o 

(A-7) 

As  for  app  in  (A-2),  it  can  be  written  as 

a„  =  aiz' )  [kl  Koiz')Y'  exp  [ }  P^(z” )  dz”  ]  /  [P_(z’ )  -  P,(z’ )]  if  z  >  z’ 

z' 

=  a{z' )  [kl  Koiz'  )]■'  exp  [  -  f  P.fz*  ’ )  dz"  ]  /  [p_(z’ )  -  p^(z’ )]  if  z  <  z* 

X 

(A-8) 

The  P  ‘s  in  (A-8)  are  governed  by  the  differential  equation 


dp*  /  dz  +  p*  -  b(z)  p*  +  a(z)  =  0 

with 


(A-9) 
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Eq.  (A-9)  is  subjected  to  the  boundary  conditions 


(A-10) 


M0)  = 


;a(0) 


Kom^ld-p^ 

= - ;  V  «(-«“) 

At  a  point  of  dielectric  constant  discontinuity  (z  =  -  d )  the  ^  ‘s  satisfy  the  condition 


(A-11) 


(X./<I)P.  Lh.  =  » 

(A-12) 

The  other  elements  of  Green’s  function  tensor  (ap*,  a*p,  Oa)  can  be  related  to  app  through 
the  relation 


aij  =  /v(z»z' )  a„i2,z'  )-[l/kl  Kc]  5(z  -  z' )  Si 

(A-13) 

The  factor  fij{z,z' )  in  (A-13)  varies  fix)m  one  Green’s  function  component  to  another. 
Also  for  the  same  Green’s  function  component,  the  value  of  this  factor  may  vary  depending  on 
the  relative  locations  of  the  field  z  and  source  z'  points  to  each  other.  Explicit  expressions  for 
apx,a,p,and  Oa  are  given  in  [1]. 


For  a  field  point  located  in  the  upper  half  space  (z  >  0) ,  the  Green’s  function  components 
in  (A-1)  and  (A-2)  can  be  written  as 
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f (z,2’  ,k.rky) s  A  (z’  ,k.,ky) cxp(;fc,z) 


(A-15) 
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APPENDIX  B 


THE  AUXILIARY  PROPAGATION  FUNCTIONS 

A  complete  description  for  the  Green’s  tensor  elements  requires  aq>licit  expressions  for 
auxiliary  propagation  functions  a’s  and  P’s  (  ±  is  dn^^ied  for  convenience).  Also  it  requires 
the  evaluation  of  an  exponential  raised  to  a  power  equal  to  the  integral  of  these  auxiliary 
functions. 

To  obtain  an  explicit  expression  for  a  we  write  the  differential  equation  (A-4)  describing 

a  as 

-  j. 

a(r)  +  a*(z) 

(B-1) 

and  assume  that  over  the  interval  of  interest,  the  dielectric  properties  do  not  vary  so  much  that  the 
dielectric  constant  may  by  replaced  by  its  mean  value.  Then  integrating  both  sides  of  (B-1)  we 
obtain 

a(z)  =  Va  tan  [Va  (c  -  z)] 

(B-2) 


where  a  is  the  mean  value  of  <z(z)  in  the  interval.  In  (B-2)  c  is  a  constant  which  can  be 
determined  by  writing  (B-2)  as 

a(z)  =  Va  tan  [Vfl  (zi  -  z)  +  Vfl  (c  -  zi )] 

(B-3) 
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If  we  let  a«ai  at  z  =  zt,wecanwritea£roin(B-3)as 


a(z)a 


Va  tan  [Va(zi -  z)]  +  ai 

1-  '^tan['\/a(zi-z)] 

Va 


To  evaluate  the  exponential  in  (A-3)  we  use  (B-2)  to  write 


/a 


dz  =  -  jVa  tan[Va(z -c)] dz  =  In 
>1 


"cos[>/a(22-c)3 

^cos[Va(zi-c)] 


(B-4) 


(B-5) 


where  In  is  the  natural  logarithm.  From  (B-5),  the  exponential  in  (A-4)  can  be  written  as 


exp(  jadz)  = 


cos  [Va(22  -■  zi)  +  Va(c  -  zi)] 
cos  bJaic  -  zi)] 


(B-6) 


With  some  mathematical  manqmlations  (B-6)  can  be  reduced  to 


exp(  Ja  dz )  =  cos[Va(za  -  zi)]  +  sin[Va  (za  -  zi)] 

r,  Va 


(B-7) 


To  obtain  e^licit  expression  for  the  auxiliary  function  P  first  we  reduce  the  differential 
equation  (A-9)  governing  P  to  a  form  similar  to  the  differential  equation  (A  -  4)  governing  a . 
This  can  be  achieved  through  scaling  p  as 


Then  substituting  (B-8)  in  (A-9)  to  get 


X  +  P'*  X*  +  a(z)  +  p’ /  dz  -  b(z)  p’ X  =  0 
In  (B-9)  letting 

dk!  dz^  biz)  X 

and  dividing  the  resultant  over  X^  the  differential  equation  governing 
dp/du+p'^  +  v^O 

Which  is  similar  to  Eq.  (A-4)  describing  a .  In  (B-1 1)  we  have, 

u-kz 

By  analogy  to  (B-4),  and  (B-7)  we  can  write 

Q,  Vv  tan [Vv  X(zi -  z)  +  p, 

P  * - K - ^ - 

l-4^tan[VV  X(r,-z)] 


(B-14) 


exp[  ^  P  dz  ]  *  cos[Vv  X(z2  -  21)]  +  sin[Vv  ^,(22  -  zi)] 

n  V  V 


To  obtain  die  explicit  expression  of  X  and  v  we  substitute  for  b  form  (A-10)  into  (B-9)  to 


get 


1 

kJ'K-p^ 


dK 


(B-15) 


Integrating  both  sides  of  (B-15)  leads  to; 

ko^K  'k„^K 


(B-16) 


Then  by  substituting  (B-16)  into  (B-12)  we  get 
kc*K^ 

a(z) 


(B-17) 
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APPENDIX  C 


ANALYTIC  EVALUATION  OF  THE  INTEGRALS  IN  EQS.  (23)  and  (35) 

In  this  Aj^jendbt,  the  integrals  and  Ipp  in  Eqs.  (23)  and  (35)  will  be  evaluated 
analytically.  Since  there  is  analogy  between  a  and  P  as  has  been  discussed  in  Appendix  B,  only 
the  evaluation  of  the  integral  will  be  considered  in  detail.  Then  by  analogy  the  integral  Ipp 

will  be  evaluated. 

To  evaluate  the  integral  I^^ ,  fiom  (B-7)  we  recall 

exp(  fadz)-  cos(Vn(z2  -  zi)]  +  sin[Vn(z2  -  zi)] 

I,  Va 

(C-1) 


Also  we  write  the  power  ^ctrum  for  the  exponential  correlation  function  (Eq.20)  as 


H^(|l.|,p)  =  <|5f>y/(|B|,p) 

(C-2) 


where 


/(|  “  1.P)  =  P'^  [1  +  P  M 1  exp(-P  H ) 
P  =  [l  +  f"p*fVf 


Then  we  introduce  the  integrals 


(C-3) 
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Zcip,a)  *  j  du  /(|«|  ,p)  cos(Va  u)  -  — - - y 

[p  +a] 

z,ip,a)  =  f  du  f{\u\  ,p)  sin  (Va  u)  =  — { 3p*  +  a } 

P  [p  +<?] 

(C-4) 


For  z  >  zi  we  can  write  as 

/♦»-  [zc(p,a)  +  ^^4^  z.ip,a)] 

£  Vn 

(C-5) 

For  z<zi  we  can  write  as 

[z,(p,a)  +  +  Jd«sin(V^«)] 

(C^) 


By  evaluating  the  integral  in  (C-6)  we  can  write  as 


/4«* 


[zc(p,n)  + 


a-^(zi) 

Va 


z,(p,<i)  + 


a+(zi)-  a_(zi) 


Zfip-a)] 


(C-7) 


where 


'  Fr/r,T  -  D  ] . 

[  sin{Vtf  (zi  “  z)}  +  cos{Va(zi  -  r))]  -  sin{Va(zi  -  z))} 


(C-8) 
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By  analogy  to  (C-5)  and  (C-7)  we  can  write  the  integral  7pp  as 


/  -  <Uf>  r  .-..X  .  P-(g») 

/pp  “■  . [  2c(.P»^/  r~.  2i(p»v)  J 


A,Vv 


For  z>z\,  and 


For  z  <  zi  where, 


r/(P.v)=p-.^^^—  «p(-p|z  -  r.| )  {[  +  p(z,  -  r)  ] 


[  ^  sin{Vv(2t  -  2)}  +  cos(Vv(zi  -  z)}]  -  ^  sin{Vv(zi  -  z)}} 


(C-9) 


(C-10) 


(C-11) 


All  parameters  in  (C-9),  (C-10)  and  (C-1 1)  are  given  in  (B-6),  (B-7),  (B-16),  and  (B-17). 
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